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Abstract 

We present a practical solution to the "sign problem" in the auxiliary field 

Monte Carlo approach to the nuclear shell model. The method is based on 

extrapolation from a continuous family of problem- free Hamiltonians. To 

demonstrate the resultant ability to treat large shell- model problems, we 

present results for 54 Fe in the full /p-shell basis using the Brown-Richter 

interaction. We find the Gamow- Teller (3 + strength to be quenched by 58% 

relative to the single-particle estimate, in better agreement with experiment 

than previous estimates based on truncated bases. 
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Recent publications have described quantum Monte Carlo methods for exact so- 
lution of the nuclear shell model. The methods are based on the Hubbard- Stratonovich 
(HS) representation || of the imaginary-time many-body propagator in terms of one-body 
propagators of non-interacting nucleons moving in a fluctuating field. Thermal averages can 
be calculated, as can ground state properties; errors arise only from discretization and sta- 
tistical sampling, both of which can be controlled. As these computations scale much more 
gently with the number of single particle orbits (N s ) and/or the number of valence nucleons 
(N v ) than do direct diagonalization techniques, they hold great promise for treating very 
large model spaces. 

Unfortunately, the applicability of shell model Monte Carlo calculations has heretofore 
been limited by the "sign problem" generic to all fermionic Monte Carlo techniques |j],[|fi||| . 
The sign of the integrand may vary from sample to sample and the net integral results 
from a delicate cancellation that is difficult to reproduce with a finite number of samples. 
The problem is well-documented (and as yet unsolved) in simulations of correlated electron 
systems [Q. Except for an important, yet schematic, class of nuclear interactions 0, we 
have found that all realistic nuclear shell model Hamiltonians suffer from a sign problem. 

In this Letter, we report a practical solution to the sign problem and present the first 
realistic calculation of a mid- /jo-shell nucleus, 54 Fe ||. Our method is based on an extrap- 
olation of observables calculated for a "nearby" family of Hamiltonians whose integrands 
have a positive sign. Success depends crucially upon the degree of extrapolation required. 
We have found that, for all of the many realistic interactions tested in the sd- and /p-shells, 
the extrapolation required is modest, amounting to a factor-of-two variation in the isovector 
monopole pairing strength. 

A general time-reversal invariant Hamiltonian with two-body interactions can be brought 
to the form 

H = J2 (e* a O a + e a O a ) {O a , O a ) , (1) 

a * a 

where the O a are a convenient set of one-body operators and O denotes the time-reverse of 
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O. For real V a , H in Eq. (1) is a manifestly time-reversal invariant. The auxiliary field Monte 
Carlo approach utilizes the HS representation of the imaginary-time many-body propagator 
U = exp(—(3H) as a path integral over one-body propagators in fluctuating auxiliary fields. 
Upon introducing N t time slices of duration Af5 = f3/N t and complex c-number auxiliary 
fields a an (n = 1, . . . , N t ), we can write the canonical expectation value of an observable O 
as 

= Tr (Oe-W) ^ jD[a}W{aMa){Q) a 
{ Tr(e-P H ) ~ J D[a]W(a)$(a) ' { ) 

Here, the approximation becomes exact as 

N t — > oo and the metric is D[a] = Tl a ^ n [da an da* an /S.f3\V a \/2ir]. The non-negative weight 
is W(a) = C( cr ) ex P( — S \Va\ Wan\ 2 Aj3), where ((a) = TrU a is the canonical partition func- 
tion of the one-body evolution operator U a = Ujy t . . . Ui, where U n = exp(—A(3h n ), and the 
one-body Hamiltonian for the n th time slice is h n = S a (e* + s a V a a an )O a + (e a + s a V a a^ n )O a , 
with s a = ±1 for V a < and s a = ±i for V a > 0. The "sign" is $(cr) = C(^)/|C(^)I and 
((D)^ = Tr (OUo) I ({&)■ Both ((a) and (C?) CT can be evaluated in terms of the N s x N s matrix 
U CT that represents the evolution operator JJ a in the single-particle space. 

The sign problem arises because the one-body partition function ((a) is not necessarily 
positive, so that the Monte Carlo uncertainty in the denominator of Eq. (2) (the ^-weighted 
average sign, (<&)) can become comparable to or larger than ($) itself. In most cases ($) 
decreases exponentially with (3 or with the number of time slices ||. 

An important class of interactions free from the sign problem (i.e., $(cr) = 1) was found 
in Ref. 0. This occurs when V a < for all a in Eq. (1). In that case, s a = 1 for all a, 
so that both h n and U a are time-reversal invariant. The eigenvectors of U ff then occur as 
time-reversed pairs with complex conjugate eigenvalues Aj, A* (i — 1, . . . , N s /2), the grand 
canonical partition function ((a) = IL|1-|-Aj| 2 is positive definite, and the canonical partition 
function for even N v is also positive definite. 

Based on the above observation, it is possible to decompose H into its "good" and "bad" 
parts, H = He + Hb, with 
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E G = J2(e* a 6 a + e a O a ) + - £ V a {O a , O a } 

a Z V a <0 

r b = \T, V a {O a ,O a ) . (3) 

1 V a >0 

The "good" Hamiltonian Hq includes, in addition to the one-body terms, all the two-body 
interactions with V a < 0, while the "bad" Hamiltonian Hb contains all interactions with 
V a > 0. By construction, calculations with Hg alone have <&(er) = 1 and are thus free of the 
sign problem. 

We define a family of Hamiltonians H g that depend on a continuous real parameter g as 
H g = Hg + qHb, so that H g= % = H. If the V a that are large in magnitude are "good", we 
expect that H g=0 = Hg is a reasonable starting point for the calculation of an observable 
(O). One might then hope to calculate (0) g = Tr (Oe _/3H9 )/Tr (e _/3 ^ 9 ) for small g > and 
then to extrapolate to g — 1, but typically (<&) collapses even for small positive g. However, 
it is evident from our construction that H g is characterized by $(c) = 1 for any g < 0, 
since all the "bad" V a (> 0) are replaced by "good" gV a < 0. We can therefore calculate 
(0) g for any g < by a Monte Carlo sampling that is free of the sign problem. If (0) g is a 
smooth function of g, it should then be possible to extrapolate to g — 1 (i.e., to the original 
Hamiltonian) from g < 0. We emphasize that g = is not expected to be a singular point 
of (C?) 9 ; it is special only in the Monte Carlo evaluation. 

In the nuclear shell model, the two-body interaction can be written in a density decom- 
position as 

9 E^ T (ac,bd)^2(-) M p K MT(ac)p K -MT(bd) . 

1 abcdKTir M 

Here p KM r = P P KM + {~) T Pkm ( t = °> !)> Pkm ( ac ) = ( a l x a c ) K M is the one-body density 
operator for the pair of proton or neutron orbits (a, c) coupled to angular momentum K 
and its ^-projection M, and ir = (— Y a+lc = ^—yb+ l d j s fa e parity. The matrices E\ T are 
constructed from the two-body matrix elements VJ T (ab,cd) of good angular momentum J, 
isospin T, and parity 7r through a Pandya transformation. For interactions that are time- 
reversal invariant and conserve parity, the E\ T {i, j) are real symmetric matrices that can be 

4 



diagonalized by a real orthogonal transformation. The eigenvectors Pkm{oc) play the role of 
O a in Eq. (1), and the eigenvalues A^ vr (a) are proportional to V a . In the Condon-Shortley 
|7|] convention Pkm = ir(—) K+M Pk-m so that the "good" eigenvalues satisfy sign [A^ 7r (a)] = 
7i(—) K+1 ||. To minimize the number of auxiliary fields required, we use the freedom to add 
an arbitrary symmetric interaction to H [[J and choose Vf T=1 = VJ^ =0 so that Ext=i = 0. 
Ekt=o is then uniquely determined by the anti-symmetric part of the interaction through 
the combination \ Vjr = o + K/r=i) • 

To demonstrate the viability and utility of the method, we have applied it to the mid- /jo 
shell nucleus 54 Fe using the realistic Brown- Richter interaction || . The number of m-scheme 
Slater determinants describing the 6 valence protons and 8 valence neutrons moving among 
the N s = 20 single-particle states of the 0/7/2,5/2 an d 0p 3 / 2) i/2 orbitals is ( 2 6 °) ( 2 g °) ~ 5 • 10 9 . 
For comparison, the largest model space treated by standard diagonalization techniques is 
currently 48 Ti [|10| where the m-scheme dimension is ~ 7 • 10 6 . 

Figure 1 (upper) shows the eigenvalues Vx n a = 7r ( — ) K ^Kn{ot) of the Brown-Richter 
interaction; only about half of the eigenvalues are negative. However, those with the largest 
magnitude are all "good" . It is possible to use an inverse Pandya transformation to calculate 
the usual two-body matrix elements Vj T (ab, cd) for the "good" and "bad" interactions, 
allowing the matrix elements of Hq to be compared in Fig. 1 (lower) with those of the full 
interaction. The greatest deviation is for J — 0, T — 1 (the monopole pairing interaction), 
where Hq is about twice as attractive as the physical H . In all other channels, Hq and H 
are quite similar. 

We have performed Monte Carlo calculations for (3 = 2 MeV -1 using N t = 32 (so that 
A(3 = 0.0625 MeV" 1 ). For g = -1,-0.8,-0.6,-0.4,-0.2, and 0, we took approximately 
3300 uncorrelated samples. The computations were performed on the Intel Touchstone 
DELTA 512-node parallel computer, where each node is an Intel i860 processor. Each node 
produced and analyzed a sample in about 4 minutes, so that each value of g took about 
25 minutes in total. Selected calculations for larger values of (3 or N t show that we have 
converged to the true ground-state properties. 

5 



The results for various observables are shown in Fig. 2. The extrapolations to the physical 
Hamiltonian (g = 1) are done by least-squares polynomials. For each observable except 
(H), the degree of the polynomial is chosen to be the lowest for which \ 2 P er degree of 
freedom is less than 1; linear or quadratic extrapolations are almost always sufficient. For 
(H), the variational principle implies the additional constraint of vanishing derivative at 
g = 1, in which case a quadratic or cubic polynomial is used. We have also calculated 
response functions R(t) = (O\t)O(0)) by polynomial extrapolation of our calculations 
of \n[R g (r) / R g (0)} for g < 0. Fitting ln[i?i(r)/i?i(0)] to a polynomial in r allows us to 
determine moments of the normalized strength function fo{E), such as E = J Efo{E)dE. 
Our overall method was checked in detail [|ll|] by comparison with direct diagonalization 



in the s<i-sheil using the Brown- Wildenthal interaction [|12j and in the lower /p-shell ( 44 Ti) 
using the Brown- Richter interaction ||. 

Table 1 summarizes the extrapolated results for various observables. Note that the 
statistical uncertainty in these values is proportional to the uncertainties in the Monte 
Carlo results for g < 0, and so can be reduced by increasing the number of samples. The 
calculated first moment of the isoscalar quadrupole strength function, (1.25 ± 0.16) MeV, 
should be compared with the empirical excitation energy of the first 2 + state, 1.408 MeV. Our 
estimate for the B(E2) for the decay of this state assuming free nucleon charges (and that 
this transition has all of the strength) is (96 ± 1) e 2 fm 4 , while effective charges (e p , e n ) = 
(l.l,0.1)e would be required to reproduce the experimental value of 126 e 2 fm 4 . These 



charges are significantly smaller than the (1.35, 0.35)e used in truncated calculations |T3 
or the (1.33, 0.64)e used in the lower fp shell ||. The total mass quadrupole strength, 
(Q 2 ) = (1482 ± 84) fm 4 , is significantly larger than the simple single particle estimate of 
380 fm 4 . The total Ml strength ((Ml) 2 ) = (14.1 ±0.4) fi 2 N is quenched relative to the single 
particle estimate of 42.55 [i 2 N . It is also interesting to note that the occupation numbers of 
the single particle orbits are smeared across the Fermi surface. 

Of particular physical interest are the Gamow- Teller operators. Our calculations exactly 
satisfy the sum rule ((GT_) 2 ) — ((GT + ) 2 ) = 3(N — Z) = 6. The single particle estimate for 
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((GT + ) 2 ) corresponding to the f 7 / 2 proton — > f 5 / 2 neutron transition is 10.28 [[H], so the shell 
model Monte Carlo value of 4.32 ±0.24 is quenched by 58%. This value is comparable to the 



experimental result of 3.1 ±0.6 [15[], but significantly smaller than previous estimates of 6.40 
or 6.70 based on truncated bases |13[]. The additional quenching on the full space correlates 
with the enhanced B(E2, 2f, — > Of), (i.e., smaller effective charge), as was surmised in [1TB 



Direct comparison with experimental Gamow- Teller strength functions requires that we 
know the energy of the daughter ground state relative to 54 Fe. Since the 54 Co ground state 
is the isobaric analog state (IAS) of the 54 Fe ground state, we find a mean (p, n) excitation 
energy of E x = (6.13 ± 0.17) MeV. This is in agreement with the systematics of Nakayama 
et al. ]TJJ, which give E GT - — = 5.81 MeV, but is somewhat low relative to the 

experimental value of 8.2 MeV [p75 ] . When our calculations of the mean (n,p) excitation 
energy are corrected for the Coulomb energy (including exchange) and the nucleon mass 
difference, we find E x = (1.24 ± 0.2) MeV, to be compared with the experimental centroid 



of 3 MeV ||15[| . A more consistent theoretical value of E x can be obtained by calculating the 
mass differences of the A = 54 isobars within the shell model Monte Carlo JTTJ . 

The method presented in this Letter is a practical solution to the sign problem for realistic 
shell model interactions. A full-basis calculation of 54 Fe with the Brown-Richter interaction 
shows the feasibility of the method, with significant quenching of the Gamow- Teller (5 + 
strength. Systematic studies of the temperature, nuclide, and interaction dependence of 
these calculations will be reported elsewhere. Our techniques also enable the determination 
of an optimal effective interaction and effective operators in a greatly enlarged model space. 
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FIGURES 

FIG. 1. Upper: The eigenvalues V a of the Brown-Richter interaction in the /p-shell. Eigen- 
values for each particle-hole angular momentum K are plotted in increasing order. Bottom: The 
two-body matrix elements VjT=i(ab, cd) of the Brown-Richter interaction (solid circles) and its 
"good" part (open circles), for J < 4. The ordering for each J is arbitary. Plots of the Vjt=o and 
the remaining T = 1 matrix elements (not shown) are similar to those shown for J > 1. 

FIG. 2. The results of the Monte Carlo calculations for 54 Fe at (3 = 2 MeV -1 for several 
observables as a function of g < 0. Q — Qp + Q n is the isoscalar quadrupole, Q v — Qp — Q n 
is the isovector quadrupole, GT + is the Gamow- Teller operator changing a proton to a neutron, 
and Ml is the magnetic moment operator using the free-nucleon g-factors. The lines are polyno- 
mial extrapolations; the extrapolated values and corresponding uncertainties are shown at g = 1. 
The extrapolation is linear for (Ml 2 ), but quadratic for (Q 2 ), (Q 2 ), and (GT 2 }. For (H), the 
extrapolation is cubic with the constraint of vanishing derivative at g = 1. 



10 



TABLES 



TABLE I. Monte Carlo results for 54 Fe. 



(H) = -55.5 ± 0.5 MeV 

Total strength E (MeV) 



Isoscalar quadrupole 




(Q 2 ) = 


1482 ± 84 fm 4 


1.25 ± 0.16 


Isovector quadrupole 




(Ql) = 


381.3 ± 33.8 fm 4 


12.7 ± 0.2 


Gamow- Teller (p, n) 




<(GT_) 2 > = 


10.32 ±0.24 


6 13 4- n 17 


Gamow- Teller (n,p) 




<(GT + ) 2 ) = 


4.32 ±0.24 


9.7 ± 0.2 


Ml 




((Ml) 2 ) = 


14.1 ±0.4 n% 


8.6 ± 0.7 






Occupation Numbers 




Protons 




Neutrons 




(a)a) h/2 


= 4.92 ± 0.03 




( a t a ) /7/2 =6.35 ±0.03 




( ata )p 3 /2 


= 0.56 ±0.02 




( a ta} P3/2 = 0.86 ± 0.02 




( at a}pi/ 2 


= 0.11 ±0.01 




(at a } Pl/2 = 0.17 ±0.01 




( ata )/ 5/2 


= 0.41 ±0.01 




(ata} /5/2 = 0.61 ±0.01 
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